Proximity effect in the presence of Coulomb interaction and magnetic field 
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We consider a small metallic grain coupled to a superconductor by a tunnel contact. We study 
the interplay between proximity and charging effects in the presence of the external magnetic field. 
Employing the adiabatic approximation we develop a self-consistent theory valid for an arbitrary 
ratio of proximity and Coulomb strength. The magnetic field suppresses the proximity-induced 
minigap in an unusual way. We find the phase diagram of the grain in the charging energy - 
magnetic field plane. Two distinct states exist with different values and magnetic field dependences 
of the minigap. The first-order phase transition occurs between these two minigapped states. The 
transition to the gapless state may occur by the first- or second-order mechanism depending on 
the charging energy. We also calculate the tunneling density of states in the grain. The energy 
dependence of this quantity demonstrates two different gaps corresponding to the Coulomb and 
proximity effects. These gaps may be separated in sufficiently high magnetic field. 

PACS numbers: 74.45.+C, 74.50.+r, 73.21.-b 



I. INTRODUCTION 

A normal metal in contact with a superconducting 
lead acquires some superconducting properties; this phe- 
nomenon is known as the proximity effect (for review see, 
e.g., Ref. 01). In particular, the electron spectrum of the 
normal metal becomes gapped; at low interface trans- 
parency this gap is much smaller than superconducting 
A, hence the name minigap. This phenomenon is due to 
the Cooper pairs of the superconductor that penetrate 
into the normal metal and induce weak superconductive 
correlations there. 

If a normal metal part of the superconductor - nor- 
mal metal junction is a small grain, then the Coulomb 
effects come into play (adding an electron to the grain 
costs charging energy). They are mostly pronounced 
in tunneling experiments when the differential conduc- 
tance between the normal external tip and the small 
grain is measured [this quantity is proportional to what 
is called the tunneling density of states (TDOS)]. The 
Coulomb repulsion between tunneling electrons reduces 
the current. This phenomenon is known as the tunnel- 
ing anomaly or, in the zero-dimensional case, Coulomb 
blockade* 2 " Except for the charging energy, an essential 
parameter governing the efficiency of the Coulomb block- 
ade is the interface conductance G between the grain and 
the lead. In the case of a normal lead the Coulomb block- 
ade is developed^ at G <C 1 (we measure G in units of 
e 2 /h) and disappears? at G 1 due to the fact that an 
electron tunneling to the grain is rapidly transferred to 
the lead, thus not blocking the tunneling of the next elec- 
tron. At the same time, the Coulomb blockade persists 
even at G 3> 1 if the lead is superconducting because 
a single electron cannot escape into the lead due to the 
gap A in its single-particle density of states (DOS). This 
situation was studied by Matveev and Glazman. 4 

Apart from the tunnel current, the Coulomb interac- 
tion suppresses the proximity effect as well. The charg- 



ing energy prevents Cooper pairs from tunneling to and 
from the normal grain and thus destroys the supercon- 
ducting order induced on the normal side of the junction. 
Qualitatively, the Coulomb interaction is trying to fix the 
charge (neutrality) of the normal grain while the proxim- 
ity effect fixes the phase. Charge and phase are conjugate 
variables obeying the uncertainty principle: they may not 
be fixed simultaneously. Recently, a full quantitative de- 
scription of this competition between the proximity and 
Coulomb effects in superconductor - normal metal junc- 
tions at large interface conductance was derived in Ref.0. 
Obviously, the Coulomb interaction diminishes the mini- 
gap. In this paper we extend the results of Ref. la, includ- 
ing the external magnetic field into consideration. The 
magnetic field can easily be varied in experiment, and we 
demonstrate that the minigap is qualitatively sensitive to 
the strength of the magnetic field. 

We consider a normal grain connected to a supercon- 
ducting lead by tunnel junctions (low interface trans- 
parency) with large conductance G (determined by the 
product of transparency by the number of channels). We 
assume the zero-dimensional limit; i.e., the Thouless en- 
ergy £/rh = D/d 2 is larger than all other relevant en- 
ergy scales of the system including the superconductive 
gap A in the leads (here D is the diffusion constant in 
the grain and d is its characteristic size). In this limit, 
the proximity-induced minigap^ is E g = GS/4, provided 
Eg <C A and S being the mean level spacing per one spin 
projection in the grain. 

Our aim is to take into account effects of Coulomb in- 
teraction and magnetic field. The characteristic Coulomb 
energy Ec — e 2 /2C is assumed (similarly to Ref. 0) to 
lie in the same range as E g : 
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The capacitance C of the grain already takes into ac- 
count the renormalization C = Cq + e 2 G/2A due to vir- 
tual quasiparticle tunneling^ This renormalization as- 
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sures the inequality Ec -C A and allows arbitrary small 
geometric capacitance Cq. 

We consider relatively weak magnetic fields H that do 
not affect the superconducting lead; i.e., H should be 
much smaller than the critical field of the lead: H <C 
H C 2 = $o/27r£|, where <!>o is the flux quantum and £s = 
\J Ds/ A is the superconductive coherence length. Note 
that since the paramagnetic limit is much larger than 
H C 2 , the condition H -C A is certainly fulfilled (here and 
below we express H in energy units dropping the factor 
9Hb/2). 

If the normal grain is sufficiently small, we can ne- 
glect the orbital effect of the magnetic field in the grain. 
Indeed, as d decreases, the critical field due to orbital ef- 
fects in a superconductor with a gap E g grows as $o/£d, 
where £ = y/D/E g is the coherence length corresponding 
to the proximity-induced superconductivity. At the same 
time, the Zeeman effect of the magnetic field determines 
the paramagnetic limit with the critical field of the order 
of Eg independent of the grain's size. Hence the disre- 
gard of the orbital effect is justified for small grains with 

d«$o/v/^- 

To take into account both the proximity and charg- 
ing effects, the adiabatic approximation was employed in 
Ref. 0. An inequality Ec ^> 5 provides the separation 
of energy scales: electronic degrees of freedom (which 
are contained in the matrix Q of the a model; see be- 
low) are "slow" compared to the characteristic frequency 
of the electric potential fluctuations. This allows one 
to calculate the renormalized (due to interaction) value 
of the minigap E g . The result of the competition be- 
tween charging and proximity effects is determined by 
a comparison of the charging energy Ec with the ef- 
fective "Josephson" energy Ej cx G 2 5\n(A/E g ). The 
latter has a clear meaning of the Josephson coupling 
energy^ between the superconductive reservoir and an 
imaginary weak superconductor with the order parame- 



ter Eg. The two limiting cases of the Coulomb versus 
proximity competition are the (i) weak Coulomb block- 
ade limit Ej Ec, when a small negative correction 
to the noninteracting minigap E g arises, and (ii) strong 
Coulomb blockade regime Ej <C Ec, when the minigap 
is exponentially suppressed. A self-consistent approach 
allowing for the magnetic field is developed in Sec. UTI 

As we show below, the magnetic field does not influ- 
ence the minigap E g if H < E g /2, while in the opposite 
case extra solutions of the model appear, resulting in a 
rich phase diagram describing different possible values of 
the minigap. An extensive study of these solutions is 
presented in Sec. IIHI 

The TDOS is also strongly affected by the magnetic 
field. This quantity is particularly interesting because it 
can be directly measured in the experiment as the dif- 
ferential conductance between the normal grain and a 
normal external tip — e.g., with the help of the tun- 
neling microscopy technique. The Coulomb interaction 
produces a drastic impact on the TDOS: as the charging 
energy increases, the proximity minigap gradually trans- 
forms into a Coulomb gap of the order of Ec- The mag- 
netic field significantly changes the energy dependence of 
the TDOS due to the spin polarization of the tunneling 
electrons. This effect is described in Sec. IIVI 



II. MODEL 

Technically, we employ the replicated zero-dimensional 
a models in imaginary time t. This model is formulated 
for calculating the disorder average of the nth power (n 
is the number of replicas) of the partition function, (Z n ). 
The standard representation of the partition function 10 
is given in terms of the coherent-state functional integral, 
Z n = f L>**D*e~ 5 [**'*], with the action 
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The fermionic two-component field $ = {ipi, con- 
sists of Grassmann (anticommuting) variables dependent 
on the space coordinate r, imaginary time r, and replica 
index a. The two-component structure of ^ corresponds 
to the Nambu-Gor'kov representation, which we need for 
studying superconductive correlations induced in the nor- 
mal grain by the proximity to the superconductor. The 
other notations used in Eq. are £ = (— iV) 2 /2m — fi, 
CAmp(r) is the potential of impurities, and tj are the Pauli 
matrices in the Nambu-Gor'kov domain. The contact to 
the lead is described by a tunnel term, which we will add 



to the action later. 

The action contains the fourth-order term due to 
the Coulomb interaction. The disorder averaging (we 
assume Gaussian 5-correlated disorder) induces another 
quartic term that mixes different replicas. These two 
terms are decoupled by the Hubbard-Stratonovich trans- 
formation with the help of a scalar variable 0" and a ma- 
trix field Qtt>- These two objects are determined in the 
space of replicas and imaginary times (or, equivalently, 
Matsubara energies). Q is also a 2 x 2 matrix in the 
Nambu-Gor'kov space. After the Hubbard-Stratonovich 
transformation, the action becomes quadratic in \P and 
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the Gaussian integration yields 
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Here e — id/dr is the Matsubara energy and Ti mp is the 
mean free time. The "Tr" symbol stands for the trace 
in all the three spaces of replicas, energies, and Nambu- 
Gor'kov, along with the integration in the real space. We 
assume the grain to be so small that the magnetic field 
has only a Zeeman but not orbital effect, d <§C <E> / \jDE g . 

The zero-dimensional approximation (d < y/D/A.) also 
implies that both Q and cj> do not vary in space and thus 
commute with £. 

In the derivation of the model J3J) we took advantage of 
the homogeneity of the magnetic field. The direction of 
the field is chosen to be the spin quantization axis. In this 
particular representation, the interaction with the field 
has only diagonal matrix elements within the Nambu- 
Gor'kov domain [see Eq. I* 1 the situation of any 

other direction of spin quantization or inhomogeneous 
magnetic field a more general model is needed with the 
twice larger Q matrix bearing also the spin indices. In 
our situation this spin-dependent Q matrix is reduced to 
a block-diagonal form describing the "up" and "down" 
spin states separately. The action © determines the full 
dynamics of one of these blocks ( "up" ) . The action for 
the second block differs from Eq. © only by the sign of 
H. Therefore, a physical quantity is given by the average 
of the two values calculated with the action ® at ±H. 
Having in mind this recipe, we derive all further results 
from the simplified model with the action |J3J. 

The a model derivation proceeds with the expansion 
of the logarithm in Eq. © in soft modes of the Q 
fieldiii These modes are concentrated at small energies, 
|e| < 1/rjmp, and lie on the manifold Q 2 = 1. At higher 
energies the matrix is diagonal, Q = f 3 sgne. Before ex- 
panding the logarithm, we have to exclude high-energy 
modes, associated with the fluctuations of the chemical 
potential, from the Q matrix. This is achieved by the 
gauge transformation- 2 . 



(4) 



with properly chosen phase K®. The matrix Q^ T i de- 
pends on the two imaginary time indices and is antiperi- 
odic on the interval [0, 1/T], The gauge transformation 
may not alter these boundary conditions; thus, we have 
to impose the restriction 



Kl jT - Kq = 2nW a , 



(5) 



with arbitrary integer W a . For a more rigorous calcu- 
lation one also has to take into account the half-integer 
values of W a . These values of the winding numbers im- 
ply that Q is periodic with respect to both imaginary 



time indices. The half-integer W a are responsible for the 
parity effect^ However, this effect is extremely weak in 
the proximity structure. To observe the parity effect, 
the two conditions should be fulfilled: (i) the minigap is 
of the order of Ec and (ii) the system is in the strong 
Coulomb blockade regime^ in which the "Coulomb stair- 
case" is well pronounced. These two requirements are 
strongly inconsistent due to the large value of the junc- 
tion's conductance G. Thus hereafter we consider only 
the integer W a . 

Assuming that Q contains only soft modes we expand 
the logarithm and obtain^ 
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Here the "Tr" operation implies trace in all the three do- 
mains of the Q field. The last term of Eq. 10 comes 
from the logarithm expansion at energies well above 
1/Timp- It corresponds to the static compressibility of 
the electron gas. Thus we have to choose K such that 
\<j) — K\ is minimized. The electric potential </>° is a real 
Bose field, 0£ = < t>Z e ~ ,WT - We separate the zeroth 

Fourier component into the integer and fractional parts, 
0« =o = 2n(W a + w a ), and choose K* to be 



C a + 2nTW a T — T — 

it, 



(7) 



The gauge transformation with such a definition of the 
phase K was proposed in Ref . [l^. Note that we still have 
the freedom of adding an arbitrary time- independent con- 
stant C a to K. Indeed, this will neither change \<j> — K\ 
nor violate the restriction JSJ. The value of this constant 
will be fixed later after we add the boundary term to the 
action. 

Now we can rewrite the action (j5J) in terms of K getting 
rid of the potential <fi. Then the functional integration 
over (j) is replaced by integration over K restricted by Eq. 
(JSJ) along with summation over W and integration over 
w in the interval [—1/2, 1/2]. We also use the inequality 
Ec ^S> S and obtain 



S[Q,K] = 



7 Tr 

(w a ) 2 



((e + iH)f 3 
W a w a ' 



2E C 



2nTw)Q 

-1/T 



E 



dr 



(K a T f 
4E C 



(8) 



Up to now we have not taken into account the tun- 
neling of electrons to and from the superconducting lead 
attached to the grain. The action |JSJ has only one non- 
trivial, but still diagonal in Nambu-Gor'kov space, term 
(the one containing f 3 ). The coupling to the supercon- 
ductor will induce off-diagonal contributions as well. To 
derive the a model with the boundary termii one has 
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to add the tunneling term to the single-particle Hamil- 
tonian. This term will appear in the argument of the 
logarithm in Eq. Expansion to the second order in 
the tunneling amplitude then leads to an additional con- 
tribution to the action, 

S b = ~Tr(Q 1 Q a ) = ~Tr{Q 



if 3 K Q e -i? 3 K 



(9) 

which is to be added to Eq. (JHJ). Generally, the bound- 
ary term contains the trace of the product of Q matrices 
on both sides of the contact (we denote them Qi and 
Q2)- Varying the full action, including the boundary 
term, with respect to Q, one can easily obtain the well- 
known Kupriyanov-Lukichev boundary conditions^ In 
the particular case of the normal grain coupled to the 
superconductor, we express Q in terms of Q and K ac- 
cording to Eq. Q). In the bulk superconductive lead the 
Q matrix takes the value 

Qs = 2Trd ab 6(e-e')T 1 . (10) 

This form of Qs is valid at low energies e <C A. Below 
we consider various properties of the normal grain at en- 
ergies not larger than E g ; thus, the above approximation 
is suitable for our purposes. The high-energy contribu- 
tion to the action (jHJl, taking into account the energy 
dependence of Qs, leads to the renormalization of the 
capacitance C — Co + e 2 G/2A as described in Ref. 0. 

Throughout the paper we assume that the temperature 
lies in the range 



S < T < E a 



(11) 



According to Eq. JHJl, the condition T > 8 fixes w a = 0. 
The fractional part of the windin g n umber freezes at 
the same temperature as in Ref. Il4l However, con- 
trary to the normal granular system with large inter- 
grain conductance, the integer part of the winding num- 
bers may still strongly fluctuate. The reason for this is 
the unconventional Coulomb blockade effect due to the 
superconductivity^ as was discussed in the Introduction. 
At temperatures below 5, the fractional parts of the wind- 
ing numbers are no longer frozen and the separation into 
W a and w a makes no sense. One should use the ap- 
proach developed in Ref. instead. We do not consider 
this limit in the present paper. 

The parameter E g , appearing in Eq. (| 1 1 11 - gives the 
characteristic energy scale of the matrix Q. The defi- 
nition of Eg will be given below [see Eq. (|17[l ]. In the 
absence of a magnetic field this parameter was found in 
Ref.H The upper bound on the temperature allows one 
to simplify further formulas, replacing all sums over the 
Matsubara energies by the corresponding integrals. Then 
the a model action takes the form 



S = --Tr 



(e + iH)%Q 



E 



dr< 



(gr) 

AE C 



nE n 



tr 



Q^(ficos2i£? + f 2 sin2i^) 



(12) 



where "tr" denotes the trace in the Nambu-Gor'kov 
space. In the absence of a magnetic field and Coulomb 
interaction, the action (|12|l is the same as for the super- 
conductive grain with the order parameter E g — G<5/4; 
therefore, E g plays the role of the bare minigap in our 
problem^ 

Due to the condition Ec S> <5, we can employ the 
adiabatic approximation. 5 Considering if as a relatively 
"fast" variable in comparison with Q, we integrate the 
action with respect to K at fixed Q. Then we come to 
the action for Q only and employ the saddle-point ap- 
proximation. The simplest saddle point of that action is 
diagonal in replicas and Matsubara energies but not in 
the Nambu-Gor'kov space. Then the condition Q = 1 
may be explicitly resolved by the parametrization 



\ab 



2t: 5 ab S{e - e') [f 3 cos 0° + n sin ( 



(13) 



The angle 9 E is the standard Usadel angleiii Generally, 
the Q matrix may also contain a fa component that is 
not present in Eq. I|13fl . This term can always be elim- 
inated by the proper choice of the constant C a in the 
definition @ of the phase K. 

All eigenvalues of Q are ±1. In each replica and at 
every Matsubara energy we have a pair of +1 and — 1 
eigenvalues. Generally, the condition Q = 1 admits an 
arbitrary distribution of the eigenvalue signs at small en- 
ergies, \e\ < l/r imp . However, the inequality T ^> 5 
allows one to neglect these unconventional saddle points. 
This is provided by the very first term of the action 112fl . 
namely, — ird^ 1 Tr(ef3<5). The minimal Matsubara en- 
ergy is 7rT; therefore, the action of those saddle points is 
larger at least by 2ir 2 T/8. The proximity effect, which 
is accounted for by other terms of the action, makes this 
estimate even stronger at low energies. We will discuss 
this issue below. 

Substituting the ansatz (|13fl into Eq. (|12f> . we find out 
that the action for K a is local in imaginary time. This 
allows us to describe the dynamics of K a by the following 
Hamiltonian: 



H a = E. 



c 



d 2 



dK 2 



2q a cos 2K 



where the parameter q a is determined by 



2E C S 



E n 



2E C S J_ 



A 



Gfesin C ". 



(14) 



(15) 



We cut off the logarithmically divergent integration at 
A since expression (|10|l . which we used for Qs, is valid 
only at £ < A. The dynamics of K is restricted by the 
condition JSJ. The summing over all integer W a results 
in the periodic boundary conditions for the eigenfunc- 
tions of the Hamiltonian 114|) on the interval [0, 2ir}. This 
Hamiltonian has an important symmetry: it commutes 
with the transformation K t—> K + w. This is related 
to the conservation of the electron parity in the grain. 
The electron number operator is h = —id/dK. In this 
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representation the first term of the Hamiltonian (|14|l is 
diagonal while the second one can change n by ±2 only. 
Physically, this property is a consequence of the Andreev 
reflection mechanism that changes the charge by ±2e. 
Another symmetry of the Hamiltonian l|14(l is the inver- 
sion K i— > —K. This is due to the particle-hole symme- 
try. It can be lifted by an external gate, which we do not 
consider in this paper. 

The adiabatic approximation relies on the facl£ that 
the characteristic frequency of the phase K fluctuations 
is much larger than E g . At the temperatures under dis- 
cussion [see Eq. Ijllf) ]. K is frozen in the ground state of 
the Hamiltonian (fH|) with the energy E (q) = Ecap{q) 1 
where ap(q) is the zeroth Mathieu characteristic value. 
The effective action for the Q matrix in terms of the an- 
gle 9 is then 



(It 



de{e + iH) cos (9^ + E (q a ) 



(16) 

The fact that the action is represented as a sum of inde- 
pendent identical contributions from each replica is due 
to the trivial in replicas ansatz . The next step is 
the saddle-point approximation. The specific form of the 
action l|16|) implies that the saddle-point value of angle 
9 is independent of the replica index. Using this fact we 
omit all replica indices hereafter. 

The variation of Eq. ifTCI) gives tan# e = E g /(e + iH), 
where the constant E g is determined by the system of 
the self-consistency equations 



Eg 
E a 



1 dE 

2E C dq : 



q = 



2A 



E c 5 Q(E g ,H) 



Here we introduce the notation 



(17) 



n(E g ,H) = max(E g ,H) + ^max 2 (E g ,H) - E 2 g . (18) 

The last equation of Eqs. (|17|) is obtained from Eq. 115|) 
where the found value of angle 9 e was substituted. 

The parameter E g has the meaning of the rcnormalized 
minigap in the thermodynamic density of states for one 
spin subband. The thermodynamic DOS itself (for spin 
up) is obtained from Q after the analytic continuation to 
the real energy E in the following way: 



Pl(E) = -Retrf 3 Q £ 



-iE+O 



(19) 



As a result, it acquires the standard BCS form shifted 
by H (while the DOS for the spin down, pi, is shifted by 
-H). The total DOS is 



P(E) = 
p BCS (E) 



l r 

2 



„BCS 



(E + H) 



„BCS 



Re cos 9 e 



-iE+O 



(E — H) 



■Re 



(20) 



\E\ 



E 2 -E\ 



(21) 



Let us now return to the discussion of unconventional 
saddle points. Suppose at some energy e and in a par- 
ticular replica the Q matrix has equal eigenvalues and 
hence is proportional to the identity matrix fo instead of 
Eq. (|13|) . This results in the effective exclusion of this 
replica-energy pair from the action ljl2|l . where only the 
combinations tr(fi j 2,3Q) are present. The parameter q, 
given by Eq. i|15|) . is also reduced by the contribution 
from the energy e. Due to Eq. i|16|) . the total change of 
the saddle-point action is 



AS 



2tt 



e cos ( 



dE nE g 
dq E C S 
_ 2tt 
~~ T 



sin 6L 



(e + iH) 2 + E 2 (22) 



The last identity is based on the self-consistency equa- 
tion H17J1 . The minigap E g is also changed; however, this 
leads to a higher-order correction in comparison with Eq. 
(I22|l . It is easy to see that the real part of AS is not less 
than 27r|e|/<5 > 2tt 2 T/S for any values of E g and H. Thus 
the estimate based on the first term of Eq. (|12[) becomes 
even stronger when the other terms are taken into ac- 
count. For example, in the case of zero magnetic field 
the lower bound is increased to 2TrE g /5. This means 
that the proximity effect gives an additional ground for 
neglecting the saddle points of the form other than Eq. 

So far, we have found the main saddle point in the 
replica-trivial sector of the a model. This result is equiv- 
alent to the direct calculation of the free energy of the 
system, averaged over disorder. Indeed, the form of the 
action lltj|) implies that the average partition function 
obeys the identity (Z n ) — (Z) n . Using this identity and 
putting the number of replicas to 1, we have for the free 
energy T = -T(lnZ) = -Tm(Z) = T S\ n=1 . Finally, 
substituting the saddle-point solution 9 e into the action 
(|l(j[) and noting that the imaginary-time integration sim- 
ply yields a 1/T multiplier in this expression, we find 



de(e + iH) 2 



(e + iH) 2 



+ Ep(q). (23) 



This free energy has the meaning of the Landau-Ginzburg 
functional while E g plays the role of the order parameter. 
The integral in the above expression contains a divergent 
contribution from high energies. As in the standard the- 
ory of superconductivity, we get rid of this divergence, 
subtracting the value of the free energy in the "normal" 
state, i.e., at E g = 0. Then the result of integration is 



N 



H 



El 



hi- 



2A 



n{E g ,H 



H - Jmvx?{E g ,H) -E 2 )+ Ep(q). (24) 
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The self-consistency equations <|17|) can be obtained by 
varying this free energy functional. The solution of Eqs. 
B17J1 gives extrema of the free energy; in particular, the 
trivial solution E g = always satisfies Eqs. <|17|) . We can 
also estimate fluctuations near the found extremal points. 
The complete calculation taking into account all possible 
fluctuating modes including those that break the replica 
symmetry is cumbersome, but leads to a simple result: 
the saddle-point approximation is valid provided E g 5. 
The details of this calculation for zero magnetic field can 
be found in Ref. Il8l where it was shown that the fluctu- 
ations produce a negligible correction to the Josephson 
current in the Coulomb-blockaded junction between two 
superconductors via the normal-metallic grain if E g S. 



III. THERMODYNAMIC MINIGAP 

Generally, the system of the self-consistency equa- 
tions Q17[) is not analytically solvable. Nevertheless, the 
ground-state energy for the Hamiltonian Ijl4(l can be ex- 
plicitly found in the two limiting cases of small and large 
q (physically, these limits correspond to the strong and 
weak Coulomb blockade, respectively). Then Eqs. i|17|) 
allow an explicit solution. 



A. Strong Coulomb blockade 

We start with the case of strong Coulomb blockade. 
This limit implies q <C 1 , which means that the Coulomb 
energy Ec is much larger than the effective Josephson en- 
ergy Ej. The phase K, which is governed by the Hamil- 
tonian (|14() . is delocalized and strongly fluctuates. At 
the same time, the charge of the grain is almost fixed. 
The Coulomb blockade wins the competition versus the 
proximity effect; the minigap is exponentially suppressed. 

At q <C 1, the potential energy in the Hamiltonian 114|) 
can be considered as a perturbation. The perturbation 
theory yields the ground-state energy E {q) = —Ecq 2 /2. 
Then the self-consistency equations (|17fl lead to 



2E C 5 



= In- 



2A 



(l(E g ,H) 



(25) 



Nonzero solutions exist only below the following value of 
the magnetic field: 



H* = 2Aexp 



2E C 5 \ 
El J 



(26) 



The quantity Q(E g ,H) is defined by Eq. l(T%|) and has 
different meanings for H below and above E g . As a re- 
sult, at H < H^, Eq. (|25|l has two nonzero solutions 



corresponding to these two cases: 



E„ 



y C Hl(2H-H*) 



at H < iJ c s , 

at iff /2 < H < H$. 



(27) 
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FIG. 1: E g (H) dependence in the limit of the strong Coulomb 
blockade. The solid line corresponds to the value of E g that 
gives the absolute minimum to the free energy. The first- 
order phase transition occurs at H = Hi, where E g vanishes 
abruptly. The dashed line shows other solutions of the self- 
consistency equations 11711 . These extra solutions exist only 
in the interval H^ < H < i/f , where the gapped and gap less 
states coexist. The curves are plotted for EcS/Eg = 4.5. 



Other parameters are G = 40 and A/E g 



150. 



Here the conditions for the two branches, H < E g (H) 
and H > E g (H), are rewritten in terms of the fixed value 
iff. The double- valued structure of the whole solution 
becomes clear in this representation. 

One can easily see that the first branch in Eq. 
which we call the gapped (S) state, corresponds to a lo- 
cal minimum of the free energy 124fl . while the second 
branch gives a maximum. The gapless (N) state E g = 
[this trivial solution of Eq. 12511 exists at any if] min- 
imizes the free energy if H > iff j 2 and maximizes it 
otherwise. The fields iff and = are the abso- 

lute instability fields for the S and N states, respectively. 
In the interval < H < iff, the two minima of the free 
energy coexist. At some value of magnetic field Hi lying 
in this interval the energies of the two states are equal. 
This is the phase equilibrium point where the first-order 
phase transition occurs. Using the free energy 124fl . we 
find this critical field: 



H 



V2' 



(28) 



The E g {H) dependence is illustrated in Fig.Q] 

The mechanism underlying the first-order phase tran- 
sition between the gapped and gapless state is exactly 
the same as in a bulk ferromagnetic superconductor. 19 
The exchange field of the ferromagnet plays the same 
role as the magnetic field in our case. The correspon- 
dence becomes complete if the superconductive pairing 
constant is taken to be A = E g /2Ec6, the Debye cutoff 
ujd is replaced by A [see Eq. l(2l)|) ]. and the order pa- 
rameter is E g . Then the critical magnetic field Hi, at 
which the first-order phase transition occurs, is simply 
the Clogston-Chandrasekhar critical fields 
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B. Weak Coulomb blockade 

Now we turn to the opposite limit of a weak Coulomb 
blockade and large q, which means that the charging en- 
ergy Eq is much smaller than the effective Josephson 
energy Ej. The cosine potential in the Hamiltonian i|14|) 
strongly localizes the phase K near and 7r values. At 
the same time, the fluctuations of charge are strong. The 
proximity effect wins against the Coulomb blockade and 
the minigap is only slightly suppressed in comparison 
with its bare value E g . 

To solve Eqs. (|f 7H . we approximate the deep minima of 
the cos 2K potential by a one-dimensional oscillator with 
the ground-state energy E (q) = —2Ec{q — ^/q)- Then, 
solving Eqs. i|17|) . we find a small correction to the bare 
minigap: 

E9 = E9 ~l^AME g ,H))- (29) 

This dependence is again qualitatively different for mag- 
netic fields above and below E g (approximately equal to 
Eg). At small magnetic field the minigap does not de- 
pend on H. The state with the field-independent mini- 
gap, coinciding with the zero-field value E g (0), is similar 
to the S state in the strong-Coulomb-blockade regime. 
At higher fields E g is logarithmically diminished. This 
state will be referred to as S'. A more accurate analy- 
sis is needed to investigate the vicinity of the H = E g 
point. It turns out that the first-order phase transition 
found in the opposite limit of strong Coulomb blockade 
persists. However, now the minigap E g , being indepen- 
dent of magnetic field at low H, experiences a very small 
steplike decrease and then gradually diminishes. This is 
the first-order transition S-S'. The E g (H) dependence is 
shown in Fig. |2 where the inset illustrates details near 
the H = E g (0) point. 

From the free energy 124(1 , we straightforwardly calcu- 
late all details of the first-order phase transition, which 
now occurs between two different gapped states, in the 
limit 5> 1. We omit this bulky calculation and give the 
results only. The magnetic field H®, at which the S state 
becomes absolutely unstable, is = E g (0). The field 
of absolute instability of the S' state we denote by i/f . 
Along with the critical field H\ they are 

B% = E g {0)- K fx, Hl = E g (0)-^-x, (30) 

x = 5 , 31 

16£2ln 3 (2A/£ g ) V ' 

where E g (0) is taken from Eq. I|29|) . The parameter x, 
which determines the scale of the phase coexistence re- 
gion, is linear in small Ec but contains also an enor- 
mously small numerical coefficient. This is the reason 
why this region is extremely small in Fig. [21 

When the first-order transition occurs, the free energy 
has two minima with identical values. What is the en- 
ergy barrier between these two minima? This barrier is 



0.8 




H/E g 

FIG. 2: E g (H) dependence in the limit of the weak Coulomb 
blockade. The inset is a close-up of the H — E g (0) point 
that is encircled in the main plot. As in Fig.0 the solid line 
shows the value of E g that gives the absolute minimum to 
the free energy, while the dashed line corresponds to extra 
solutions of the self-consistency equations 1171 . The curves 
are plotted for Ec5/E g — 1.5. Other parameters are G = 40 
and A/E g = 150. 

also numerically very small and equals 2E g x 3 ^ 2 /3S. Ob- 
viously, when the height of this barrier becomes compa- 
rable with the temperature, fluctuations smear the first- 
order transition. Then a crossover between S and S' 
states occurs instead of a phase transition. 

Finally, when the magnetic field is high enough (be- 
yond the scope of our model), it suppresses the supercon- 
ductivity in the lead. In the absence of Coulomb effects, 
the minigap persists as long as the lead is superconduct- 
ing and disappears at the critical field H C 2 of the lead. 
If the weak Coulomb blockade is realized at H = 0, then 
the minigap will vanish at a field slightly smaller than 

H C 2- 

C. Intermediate case 

In this section we consider an intermediate regime 
when the Coulomb interaction is comparable with the 
proximity coupling. In Fig. [3] we present the phase dia- 
gram in the Ec~H plane. This diagram covers all limit- 
ing regimes considered above along with the intermediate 
region. 

We have already studied the first-order transition from 
S to N and S' states at large and small Ec, respectively. 
At the same time, the line of absolute instability of the S 
phase can be extracted from results of Ref. obtained at 
H = 0. Indeed, at low magnetic field H < E g , the self- 
consistency equations (|17fl do not contain H . Hence the 
minigap is independent of the magnetic field (S state) and 
coincides with the zero- field value E g (0). The maximal 
possible magnetic field for this state is — E g (0). Ob- 
viously, this result holds for any value of q. The E g (Ec) 
dependence in the absence of magnetic field was studied 
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FIG. 3: Phase diagram Ec—H. The solid line marks the first- 
order transition from the S to either N or S' state. The region 
of phase coexistence is shaded. The dashed line shows the 
second-order S'-N transition. All the three phases equilibrate 
at the triple point T. A detailed diagram in the vicinity of 
the triple point is shown in Fig. 01 The diagram is plotted for 
G = 40 and A/E g = 150. 

in Ref.H 

Now we concentrate on the S'-N second-order transi- 
tion and the vicinity of the triple point where all three 
phases equilibrate. The S'-N transition line can be cal- 
culated analytically. Any solution of the self-consistency 
equations (|17[) gives an extremum of the free energy: 
dT jdEg = 0. The normal state [E g = 0) always sat- 
isfies this condition. The normal state is stable provided 
d 2 T jdE 2 > 0. Calculating the second derivative of the 
free energy (|24[1 and then taking the limit E g — > 0, we 
easily find the critical magnetic field 

ff c N = A exp (-2E C 6/E 2 g ) . (32) 

This critical field determines the boundary of the normal 
region in the phase diagram in Fig. [3] What happens 
just below this boundary? The second derivative of the 
free energy becomes negative. If, at the same time, the 
fourth derivative is positive, then the free energy achieves 
a minimum at small E g . This is the second-order phase 
transition from the N to S' state. Otherwise, if the fourth 
derivative is also negative, then below a minimum 
at small E g vanishes and the only stable state has finite 
value of the minigap E g (Q). Thus is the normal-state 
absolute instability field for the N-S first-order transi- 
tion. In the strong-Coulomb-interaction limit, the criti- 
cal fiel d coincides with =H%/2, which we found 
in Sec. Imll 

The point on the critical line l)32|l . where the fourth 
derivative of the free energy changes its sign, is denoted 
as B (see Fig. QJ. To find this point, one should use 
a more precise value of the ground-state energy of the 
Hamiltonian 1)14(1 . taking into account the fourth-order 
perturbation correction: E (q) — Ec{—q 2 /2 + 7q 4 /128). 
Then taking the fourth derivative of the free energy (12411 , 
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FIG. 4: Phase diagram near the triple point. The solid line 
shows the first-order transition from S to N for Ec above 
the triple point and from S to S' otherwise. The dashed line 
corresponds to the critical field Hf, given by Eq. The 
dotted line denotes the absolute instability of the S' phase 
at the field Jff . This line ends at the point marked as B. 
The dash-dotted line is H^. It intersects with at the A 
point. The region of a possible metastable state is shaded. 
The diagram is plotted for G = 40 and A/E g = 150. 



we find the equation Hg = E g /7Ec<S, which, together 
with Eq. (|32|l . determines the position of the B point. 

In Fig.QJa close-up of the triple-point region is shown. 
The curve is shown by the dashed line. The line of 
absolute instability of the S' phase H~ (the dotted line 
in Fig. 0| ends in the B point. Indeed, the S' phase with 
arbitrary small E g exists only if the fourth derivative 
of the free energy is positive. Another feature of the 
phase diagram, the A point, is the point where the N- 
S' second-order transition and the absolute instability of 
the S phase occur simultaneously. Finally, between the A 
and B points on the line the triple point T lies. This 
is the point where the first- and second-order transition 
lines intersect. All three phases have the same energy in 
this point. 

Another illustration of the complicated phases struc- 
ture near the triple point is given by Fig.[5l where several 
Eg(H) dependences are shown. The leftmost curve cor- 
responds to Ec above the B point. Qualitatively, this 
case is similar to the strong-Coulomb-blockade limit (see 
Fig. ^) . The next curve is plotted for Ec taken at the B 
point. It looks much the same, but the unstable solution 
(dashed line) vanishes as the fourth rather than square 
root of the magnetic field. The next curve is for the T 
point. The first- and second-order transitions (solid and 
dashed lines) occur at the same magnetic field. The right 
but one curve is for the A point [H® — E g (0)] . The right- 
most curve illustrates the case of Ec below the A point. 
It is similar to the weak-interaction limit (see Fig. 0). 
The minigap vanishes continuously at . As Ec fur- 
ther decreases, this critical field grows exponentially and 
rapidly goes beyond the scope of our model. 
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up), we should analytically continue the expression 



FIG. 5: Minigap as a function of magnetic field (solid lines) for 
several values of Ec near the triple point. The dashed lines 
show extra solutions of the self-consistency equations 1171 . 
The curves from left to right correspond to EcS/Eg — 3.45, 
3.289 (B point), 3.152 (T point), 3.054 (A point), and 2.95. 
The other parameters are G — 40 and A/E g = 150. 



IV. TUNNELING DENSITY OF STATES 

Measuring the thermodynamic density of states is ex- 
perimentally complicated due to small size of the sam- 
ple. The tunneling technique is more practical in this 
case. The actual measured quantity is the tunnel cur- 
rent which depends on the voltage applied between the 
system and a normal-metallic external tip. The differen- 
tial conductance dl/dV extracted from this experiment 
is proportional to the local tunneling density of states 
at energy eV. The latter is determined by the imagi- 
nary part of the one-particle Green function. Without 
an interaction, the thermodynamic and tunneling DOS 
coincide. However, in an interacting system the Green 
function is "dressed" by the interaction that yields the 
difference between the two quantities. In cr-model lan- 
guage, the thermodynamic DOS is determined by the Q 
matrix [see Eq. I|19|) ]. while the tunneling DOS is given 
by a similar expression with the "dressed" matrix Q. 

In fact, the differential conductance gives the tunneling 
density of states only at zero temperature. If T > 0, the 
tunneling electrons are dispersed in the energy range of 
the order of T. The expression for dl/dV then reads 



dl 
~dV 



2R T 



dE- 



l (E + eV) 



ATcosh z {E/2T) 



(33) 



with Rt being the tunnel resistance between the tip and 
the grain. In order to measure the subtle structure of the 
tunneling DOS due to the Zeeman splitting, the temper- 
ature should be low enough: 

T <€min(E g ,H). (34) 

To calculate the TDOS for one spin projection (spin 



p t T un ( £ ) = -tr(f 3 g £e ) 

dre l£T tr ( f 3 e 4f3 ^ Q TO e~ jf3 *"° ) (35) 



to the real energies, ie — > E + iO, and take its real 
part. The angular brackets denote averaging with weight 
e -S[Q,K] w j iere t ne action is given by Eq. I|12|) . The 
TDOS for the spin down is then obtained after invert- 
ing the sign of the magnetic field. The adiabatic and 
saddle-point approximations allow one to substitute Q of 
the form l(T3|) and average over K. This procedure leads 
to the expression 



pf n ( £ ) 



duj (lo + iH ) 
» 2n < /(w+iff) 2 



-.C(e- 



UJ 



(36) 



where C(u>) is the phase correlator containing average 
over the ground state of the Hamiltonian (|14|) : 



C(u) 



dTe luJT (cos(K T - K )) 



(37) 



This quantity can be expressed in terms of the eigen- 
values E n and the eigenfunctions |n) of the Hamilto- 
nian {T1J: 



Pn 



2-Ar, 



(0|cosX|n)| 2 + |(0|sinif|n)| 2 



A n — E n — Eq, (38) 
(39) 



Now we can perform the analytic continuation to real 
energies in Eq. H36|) . In Fig. we plot the integration 
contour in the complex plain of the to variable and show 
the two poles of C(e — lj) along with the branch cut 
for the rest of the integrand. For simplicity we retain 
only a pair of poles corresponding to the nth term of the 
sum Q38|l. The analytic continuation moves the poles; as 
a result, all the singularities of the integrand reside on 
the imaginary axis. Therefore the integral along the real 
axis becomes purely imaginary and does not contribute 
to the TDOS. The real part — and hence the TDOS — is 
nonzero if a pole traverses the real axis while moving (see 
Fig- EJ) • Then the residue in this pole will determine the 
result. Note that the value of this residue is nothing but 
the thermodynamic DOS at the energy corresponding to 
the final position of the pole. To obtain the complete 
expression for the TDOS, we sum the contributions from 
all the terms in the sum l|38l) and symmetrize the result 
with respect to the spin direction: 

p t ™(E) = l[pf n (E)+p t r(E) 

= Y J PnV{\E\-A n )p(\E\-A n ). (40) 
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FIG. 6: Analytic continuation of the integrand of Eq. 136H 
in the complex ui plane. The integral is taken along the real 
axis denoted by the line with an arrow. The dashed lines are 
branch-cut discontinuities. The small dots show the poles of 
the nth term in the sum lEU for C(e-w). These poles move 
to new positions at the imaginary axis due to the analytic 
continuation. The residue in the pole that traversed the real 
axis provides the real part of the integral and determines the 
TDOS. (a) Left panel corresponds to the S phase, H < E g . 
The Coulomb and proximity gaps add up around E — 0. 
(b) Right panel corresponds to the S' phase, H > E g , where 
the two gaps are separated. 



Here -d(x) is the Heaviside step function and p(E) is the 
thermodynamic DOS given by Eq. Q2U[I. 

The resulting expression for the TDOS has a clear 
physical meaning. Every term of the sum l|40|l is obtained 
from p{E) by inserting the 2A n = 2(E n — Eq) Coulomb 
gap around E = and multiplying by the factor P n . The 
energy dependence of the TDOS is qualitatively different 
for the two gapped phases S and S', described in the pre- 
vious section. In the S phase, when H < E g , the peaks 
of the TDOS (see Fig. EJ are split due to the Zeeman ef- 
fect, but the whole picture resembles the result of Ref. 0. 
At higher magnetic field H > E g , the grain is in the S' 
phase. The Coulomb and proximity gaps are now sep- 
arated as shown in Fig. [8] The Coulomb gap is always 
centered around E = while the minigap is shifted by 
H . The density of states inside this shifted minigap is no 
longer zero due to the contribution from electrons with 
the opposite spin. 

The resulting expression l|40|l contains the matrix ele- 
ments P n and the energy level separations A n = E n — Eq . 
Both quantities can be found analytically in the lim- 
its of weak and strong Coulomb interaction. Previ- 
ously, we have mentioned two symmetries of the Hamil- 
tonian (|14fl : it commutes with operations K i— > —K 
and K >— > K + ir. Therefore, the P n coefficients are 
nonzero only for n — 4k + 1 and n = 4fc + 2. Another 
of their properties is the normalization ^2 n P n = 1, thus 
p(E) = p tun (E) = 2/5 if the energy E is far from the 
Fermi energy. The sequence P n rapidly decreases at any 
strength of the Coulomb interaction. This allows us to 
keep only n = 1 and n — 2 terms in the sum (|40|) . In the 
weak- Coulomb-blockade limit, g > 1, we approximate 
the cos 2K potential by two deep parabolic wells. The 
splitting of the two lowest levels is exponentially small. 
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FIG. 7: Energy dependence of the tunneling DOS for the S 
phase. The parameters are Ec3/Eg — 2.5, H/E g — 0.4, and 
= 150. The inset shows the details of the Zeeman 



A/E g 

splitting of the first peak, 
same manner. 
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Energy dependence of the tunneling DOS for the 
S' phase. The parameters are EcS/E g = 1.0, H/E g = 3.0, 
and A/E g — 150. The inset shows the details of the first 
peak structure. The Coulomb and proximity gap are clearly 
separated for H > E g . 



To the main order in 1/ y/q, we have 
1 



Pi = 1- 



Pi = 



1 



A l =0, A 2 = 4E C ^. 



(41) 
(42) 



In the opposite case of a strong Coulomb interaction q <C 
1, we employ the perturbation theory in q and obtain 



Pi,2=\(l±~), A h2 = E c (lTq)- 



(43) 



In the limit of a strong Coulomb interaction, at high 
magnetic fields H > H\, the minigap is absent and the 
grain is in the N phase. In this case, the TDOS ex- 
hibits the pure Coulomb gap at energies below Ec and 
is constant above this gap: p tun (E) = (2/5)tf(\E\ - E c )- 
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In fact, even i n the N regime a very small minigap of 
the order of ^S/HE g (0) [where E g (Q), the minigap at 
H = 0, coincides with iff given by Eq. J2EJ)] persists in 
the TDOS due to the fluctuations, which we neglected 
in this paper. The mechanism of this effect is similar 
to that for a superconductive grain in a strong magnetic 
field described in Ref. 0. 

V. CONCLUSIONS 

In conclusion, we developed a self-consistent theory for 
proximity and charging effects in the presence of an ex- 
ternal magnetic field. The minigap induced in the grain 
shows a complicated dependence on magnetic field. Two 
distinct minigapped states exist, and a first-order phase 
transition occurs between them. The transition to the 
gapless state is of first order from the S state and of sec- 
ond order from S'. The tunneling DOS is also different 
in S and S' states. In high magnetic field H > E g , the 
TDOS acquires two distinct gaps: a Coulomb gap at zero 
energy and a proximity minigap shifted from zero by a 
magnetic field. 

The systems discussed in this paper can be experimen- 



tally realized with the following parameters: the grain 
can be made of a noble metal and have the form of 
a disk with thickness d ~ 25 nm and diameter an or- 
der of magnitude larger. At the interface transparency 
of order 2 x 10~ 6 (insulating oxide barrier) we estimate 
S ~ 5 x 1CT 4 K, E g ~ 0.01 K, E c ~ 0.1 K, Ej ~ 1 K, and 
G ~ 100. The lead can be made of Nb; then A - 12 K 
and the magnetic fields should be smaller than 1 T. 
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